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LIMITATIONS  ON  THE  APPLICABILITY  OF  HIGH-EXPLOSIVE  CHARGES 
FOR  SIMULATING  NUCLEAR  AIRBLAST 

D,  Book.  D.  Fyfe,  M.  Picone,  Naval  Research  Laboratory, 

Washington,  D.C. 

M.  Fry,  Science  Applications,  Inc.,  McLean,  Virginia 


The  flow  fields  that  result  from  nuclear  and  high 
explosive  (HE)  detonations  are  qualitatively  alike  but 
quantitively  different.  Consequently,  care  must  be  exercised 
in  carrying  over  conclusions  drawn  from  measurements  of  HE 
tests  to  nuclear  explosions.  The  usefulness  of  HE  explosions 
for  simulating  nuclear  airblast  is  predicated  on  the  fact 
that  after  reaching  5-6  times  the  initial  radius  the  flow 
field  looks  like  that  produced  by  a  point  source  and  produces 
shock  overpressures  similar  to  those  in  the  nuclear  case. 
Numerical  simulations  of  airblast  phenomena  have  been  carried 
out  using  one-  and  two-fluid  Flux-Corrected  Transport  hydro¬ 
codes  in  one  and  two  dimensions.  The  principal  differences 
in  the  free-field  solutions  are  the  presence  in  the  HE  case 
of  a  contact  discontinuity  between  air  and  HE  products  and  of 
a  backward-facing  shock  behind  it.  Temperatures  in  the 
nuclear  fireball  are  initially  three  orders  of  magnitude 
higher;  correspondingly,  the  density  rainimum  at  the  center  of 
the  fireball  is  much  broader  and  deeper.  When  the  blast  wave 
in  a  nuclear  height-of-burst  (HOB)  situation,  undergoes 
regular  reflection  from  the  ground  only  one  peak  develops  in 
the  overpressure,  and  the  reflected  wave  propagates  upward 
rapidly  through  the  hot  underdense  fireball.  In  the  HE  case 
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part  of  the  upward-moving  reflected  wave  is  reflected 
downward  at  the  contact  surface,  producing  a  second  pressure 
peak  on  the  ground,  while  the  shock  transmitted  through  the 
contact  surface  propagates  slowly  upward.  After  transition 
to  Mach  reflection  other  differences  appear.  At  late  times 
following  shock  breakaway  the  nuclear  fireball,  unlike  the  HE 
fireball,  appears  to  develop  a  Rayleigh-Taylor  instability 
along  its  lower  edge  below  the  HOB.  The  vortices  (both 

forward  and  reverse)  are  stronger  and  form  earlier.  This  has 
important  consequences  for  fireball  rise  and  for  dust 

entrainment  and  transport  to  high  altitudes. 

Section  1 
INTRODUCTION 

In  this  paper  we  describe  a  series  of  calculations 
carried  out  as  part  of  the  ongoing  NRL  effort  aimed  at  study¬ 
ing  airblast  effects.  The  phenomena  of  chief  interest  to  us 
include  the  following:  peak  overpressures  and  pressure 

histories  on  the  ground  as  functions  of  yield,  range,  and 
height  of  burst  (HOB),  both  at  early  times  (prior  to  and 

during  transition  to  Mach  reflection)  and  at  late  times 

(after  shock  breakaway,  with  peak  pressure  in  the  range  of 
tens  of  psi);  velocity  fields,  particularly  those  associated 
with  the  toruses  (both  forward  and  reverse)  in  the  neighbor¬ 
hood  of  the  rising  fireball;  and  the  distribution  of  dust 
lifted  off  the  ground  by  the  winds  and  the  structure  of  the 
cloud  at  the  time  of  stabilization.  We  are  interested  in 
comparing  the  nuclear  and  HE  cases,  and  learning  how  much 
they  differ  from  one  another.  Our  motivation  is  to  determine 
the  extent  to  which  HE  tests  can  simulate  events  in  a  nuclear 
height-of-burst  situation. 


The  technique  we  have  employed  for  this  purpose  is 
numerical  modeling.  One-  and  two-fluid  hydrocodes  based  on 
the  Flux-Corrected  Transport  ( FCT)  shock-capturing  techni¬ 
ques1  have  been  used  to  simulate  airblast  phenomena  in  one 
and  two  dimensions.  FCT  refers  to  a  class  of  state-of-the- 
art  fluid  computational  algorithms  developed  at  NRL  in  the 
course  of  the  past  ten  years  with  supersonic  gas-dynamic 
applications  expressly  in  mind.  Simply  put,  our  procedure  is 
to  model  a  one-kton  nuclear  burst  and  its  600-ton  chemical 
equivalent,  both  at  a  HOB  of  50  m,  and  compare  the  results. 
In  order  to  validate,  initialize,  and  interpret  these  2D 
simulations,  a  number  of  ancillary  calculations  (mostly  ID) 
were  undertaken.  The  results  are  most  conveniently  exhibited 
in  terms  of  plots  of  peak  overpressure  vs  range  and  time, 
station  histories,  contour  plots  of  combustion  product  and 
total  density,  velocity  vector  plots,  and  tracer  particle 
trajectories.  Examples  of  these  are  presented  to  illustrate 
our  results  and  conclusions. 

'i.ne  plan  of  the  paper  is  as  follows:  In  the  next 
section  we  discuss  our  numerical  techniques  and  validation 
procedu^^s.  In  Section  3  we  discuss  the  free-field  (ID) 
solution  and  indicate  the  salient  differences  between  nuclear 
and  HE  cases.  Section  4  describes  the  2D  HOB  calculations 
done  for  the  HE  and  nuclear  cases.  In  Section  5  we  summarize 
our  conclusions  and  discuss  their  domain  of  validity.  We 
find  that  simulation  of  nuclear  explosions  by  HE  has  distinct 
limitations,  particularly  at  early  times,  in  the  fireball  and 
transition  regions,  and  in  the  details  of  the  dust  scouring 
process. 


Section  2 

NUMERICAL  TREATMENT 


FCT  is  a  finite-difference  technique  for  solving 
the  fluid  equations  in  problems  where  sharp  discontinuities 
arise  (e.g.,  shocks,  slip  surfaces  and  contact  surfaces).1 
It  modifies  the  linear  properties  of  a  second-  (or  higher-) 
order  algorithm  by  adding  a  diffusion  term  during  convective 
transport,  and  then  subtracting  it  out  "almost  everywhere"  in 
tne  antidiffusion  phase  of  each  time  step.  The  residual 
diffusion  is  just  large  enough  to  prevent  dispersive  ripples 
from  arising  at  the  discontinuity,  thus  ensuring  that  all 
conserved  quantities  remain  positive.  FCT  captures  shocks 
accurately  over  a  wide  range  of  parameters.  No  information 
about  the  number  cr  nature  of  the  surfaces  of  discontinuity 
need  be  provided  prior  to  initiating  the  calculation. 


The  FCT  routine  used  in  the  present  calculations, 
called  JPBFCT  (an  advanced  version  of  ETBFCT) 2 ,  consists  of  a 
flexible,  general  transport  module  which  solves  1-D  fluid 
equations  in  Cartesian,  cylindrical,  or  spherical  geometry. 
It  provides  a  finite-difference  approximation  to  conservation 
laws  in  the  general  form: 


/  $dV  »  -  /  4>(u-u  )«dA  +  /  xdA,  (1) 

6V(t)  6A( tj  6A(t 


where  <fr  represents  the  mass,  momentum,  energy  or  mass  species 
in  cell  6V(t),  u  and  Ug  represent  the  fluid  and  grid  velo¬ 
cities,  respectively,  and  x  represents  the  pressure/work 
terms.  This  formulation  allows  the  grid  to  slide  with 
respect  to  the  fluid  without  introducing  any  additional 


numerical  diffusion.  Thus,  knowing  where  the  features  of 
greatest  interest  are  located,  one  can  concentrate  fine  zones 
where  they  will  resolve  these  features  most  effectively  as 
the  system  evolves. 

The  same  transport  routine  is  employed  in  the  2D 
r-z  code  (called  FAST2D)  via  coordinate  splitting.  A  Jones- 
Wilkins-Lee  (JWL)  equation  of  state  (EOS)  was  used  for  the 
detonation  products  and  a  real -air  EOS  was  used  outside  the 
HE-air  interface.3  The  routine  was  written  in  the  form  of  a 
table  lookup,  using  interpolation  with  logarithms  to  the  base 
16  computed  by  means  of  logical  shifts.4  By  thus  taking 
account  of  the  architecture  of  the  machine  (in  these  calcula¬ 
tions,  a  32-bit-word  two-pipe  Texas  Instruments  ASC)  it  was 
possible  to  generate  very  efficient  vector  code,  decreasing 
the  time  required  for  EOS  calculations  to  a  small  fraction  of 
that  required  for  the  hydro.  The  EOS  specifies  pressure  as  a 
function  of  density  and  internal  energy.  In  mixed  cells  the 
combined  pressure  was  calculated  according  to  Dalton's  law. 

For  the  HE  calculations  the  initial  conditions  were 
taken  to  be  the  self-similar  flow  field  corresponding  to  a 
spherical  Chapman-Jouguet  detonation  at  the  time  the  detona¬ 
tion  wave  reaches  the  charge  radius  (Fig.  I)3.  This  was 
propagated  with  the  ID  spherical  code  until  the  detonation 
front  attained  a  radius  just  smaller  than  the  HOB,  at  which 
time  the  solution  was  laid  down  on  the  2D  mesh  (Fig.  2).  The 
nuclear  calculation  was  initialized  with  the  1-kton  standard5 
with  the  same  initial  radius. 

The  boundary  conditions  were  chosen  to  enforce 
perfect  reflection  on  the  ground  and  on  the  axis  of  symmetry 
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[ (<H/<2n ) Dc  “  °»  where  $  =  p,p,vfc,  and  v£c  =  0],  where  "t"  and 
"n"  denote  tangential  and  normal  components,  respectively, 
and  outflow  on  the  outer  and  the  top  boundaries  [  (d$/dn  = 
0,  where  <j>  =  p,p,vfc,vn]. 

For  the  2D  calculations  the  mesh  was  typically 
-100  x  100.  Fixed  gridding  was  used  to  minimize  numerical 
errors.  These  zones  were  2.1  m  x  2.1  m.  For  the  late  time 
calculations,  a  fixed  mesh  with  100  zones  in  the  radial  and 
200  zones  in  the  vertical  direction  was  used,  with  all  cells 
of  dimension  4.2  m  x  4.2  m. 


Section  3 

FREE-FIELD  SOLUTION 

The  well-known  Sedov  similarity  solution6  for  a 
point  blast  consists  of  a  strong  shock  (post-shock  pressure 
much  larger  than  ambient  pressure)  followed  by  a  rarefaction 
wave  (Fig.  3).  The  density  distribution  is  extremely  concave 
and  approaches  zero  at  the  origin.  The  pressure  approaches  a 
constant  as  r  •+•  0,  so  that  the  temperature  diverges  strongly. 
The  profiles  of  the  1-kton  standard  solution  (Fig.  4)  are 
qualitatively  similar,  the  temperature  being  essentially  flat 
in  rhe  fireball  region. 

The  solution  used  to  initialize  the  HE  problem, 
however,  contains  a  number  of  features  which  are  absent  in 
the  other  solutions.  These  are  most  conspicuous  in  the 
density  profile  (Fig.  2),  which  exhibits  a  contact  discon¬ 
tinuity  between  the  HE  products  and  shocked  air,  a  secondary 


distribution  behind 


shock  facing  inward  within  the  detonation  products,  and  a 
gentle  maximum  near  the  origin.  Because  of  this  last 
feature,  temperatures  are  about  three  orders  of  magnitude 
smaller  in  the  fireball  than  in  the  nuclear  case,  and  are 
nowhere  divergent. 

It  follows  that  the  speed  of  sound  in  the  nuclear 
fireball  is  much  greater  than  in  the  HE  fireball.  This  has 
two  immediate  consequences,  one  physical  and  one  numerical. 
The  first  is  that  shocks  propagating  through  the  nuclear 
fireball  travel  much  faster.  The  second  is  that  the  upper 
limit  on  the  computational  timestep,  set  by  the  Courant 
criterion 


max  6t)  <  1,  (2) 

which  usually  is  determined  by  conditions  in  the  fireball,  is 

much  smaller  relative  to  the  shock  time  scale  x  = 

s 

HOB/vshock  in  the  nuclear  case  than  in  the  HE  case. 

As  a  result,  even  though  the  leading  shock  is  well 
resolved  (over  ~  2  zones)  in  the  free-field  solution,  the 
process  of  reflection  even  at  ground  zero  (where  the  shock  is 
incident  normally)  takes  hundreds  of  timesteps.  Coupled  with 
the  property  of  FCT  (known  as  "clipping")  which  makes  the 
points  of  all  sharply-peaked  profiles  tend  to  flatten  out 
until  they  are  >  3  zones  across,  this  makes  the  rise  time  of 
the  reflected  shock  quantities  ~  10  times  longer  than  that  of 
the  incident  shock.  This  can  be  seen  using  a  ID  spherical 
model  calculation  (Fig.  4),  as  well  as  in  oblique  reflection 
in  2D.  This  spreading  is  a  problem  only  while  the  shock  is 
in  the  immediate  vicinity  of  the  reflecting  boundary..  After 
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the  reflected  shock  has  propagated  a  few  zones  back  into  the 
interior  of  the  mesh,  the  profiles  steepen  up  and  assume 
their  correct  forms  (this  has  been  shown  by  rerunning  the 
calculation  with  a  refined  mesh  and  comparing  the  peak  values 
after  reflection  with  those  predicted  theoretically). 

When  the  reflected  shock  begins  propagating  back  to 
the  origin  it  encounters  drastically  different  conditions  in 
the  HE  and  nuclear  cases.  In  the  former,  it  strikes  the 
contact  discontinuity  where  it  is  partly  transmitted  and 
partly  reflected.  The  reflected  shock  then  proceeds  outward 
until  it  reaches  the  ground  (the  end  of  the  grid  in  the  ID 
calculation),  producing  the  second  peak  in  the  station 
history  shown  in  Fig.  5.  In  contrast,  the  shock  wave 
reflected  from  the  ground  passes  unhindered  through  the 
fireball  at  high  speed  until  it  reaches  the  upper  boundary  of 
the  f-tpball  whereupon  it  reflects  back. 

As  we  shall  show,  it  is  primarily  through  these 
reverberating  shock  waves  and  the  wind  pattern  they  set  up 
that  the  HE  and  nuclear  HOB  airblasts  differ. 


Section  4 

2D  SIMUJATION  OF  AIRBLAST 


The  yield  and  HOB  (600  tons  and  166  ft,  respec¬ 
tively)  were  chosen  to  equal  the  values  used  in  the  Direct 
Course  experiment,  which  we  are  simulating.  The  Chapman- 
Jouguet  parameters  used  to  initialize  the  spherical  free- 
field  calculation  were  taken  to  be  those  for  the  NH4N03-fuel 
oil  (ANFO)  mixture  used  as  the  explosive. 


Figures  6(a)— (c)  show  the  contours  of  HE  density 
and  internal  energy  per  unit  mass  and  the  velocity  arrow  plot 
at  t=0,  just,  before  the  reflection  at  ground  zero  occurs. 
Figures  6(d)-(f)  show  the  corresponding  plots  54  ms  later, 
while  Figs.  6(g)— (i)  show  them  after  245  ms.  Note  the 
reflected  shock  proceeding  upward,  reflecting  again  off  the 
fireball,  and  propagating  back  in  a  downward  and  outward 
direction.  The  interaction  of  this  shock  with  the  radially 
inward  flow  near  the  ground  generates  the  reverse  vortex, 
which  is  clearly  seen  in  Fig.  6(i).  Note  also  the  positive 
vortex  forming  near  the  top  of  the  grid  in  the  same  plot. 
The  latter  results  when  the  upward -propagating  reflected 
shock  interacts  with  the  radially  outward  flow  near  the  top 
of  the  fireball;  it  is  not  produced  by  the  buoyant  rise  of 
the  fireball,  which  at  these  early  times  has  scarcely  begun. 

To  look  at  the  evolution  of  the  fireball  at  late 
times,  we  reinitialized  on  a  larger,  coarser  grid,  represent¬ 
ing  a  cylinder  400  m  in  radius  and  800  m  high.  The  first  300 
cycles  approximately  reproduce  the  earlv-time  results.  The 
spherical  shock  breaks  away  and  leaves  the  mesh.  The  flows 
remaining  on  the  grid  are  now  subsonic  everywhere.  Then  the 
fireball  begins  to  rise  and  the  subsequent  development  is  due 
to  the  combination  of  buoyant  rise  and  the  action  of  the 
vortices  set  up  by  the  early  shocks. 

Figures  7(a)-(b)  show  the  reaction  product  density 
and  velocities  at  0.93  s.  Note  the  "toe"  reaching  out  along 
the  ground  and  the  bulge  near  the  bottom  of  the  HE  product 
density  produced  by  the  constructive  interference  of  forward 
and  reverse  vortices.  These  features  are  accentuated  with 
the  passage  of  time;  in  Figs.  7{c)-(d)  (t  =  2.70  s),  they  are 
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even  clearer.  The  cloud  has  become  quite  elongated  verti¬ 
cally  and  shows  a  distinct  mushroom  shape.  Development  slows 
as  the  fireball  cools  and  velocities  diminish.  By  7.34  s 
(after  2600  timesteps)  the  cloud  is  almost  at  600  m.  Figures 
7(e)— (f)  show  its  form  at  this  time.  Note  that  the  maximum 
velocity  is  now  145  m/s. 

Figure  8a,  which  shows  the  trajectories  of 
passively  advected  trace*  particles  over  the  time  interval 
1.8  sec  to  3.97  sec,  displays  the  vortices  very  clearly. 
Figure  8b  shows  the  particle  paths  for  the  time  interval  3.97 
sec  to  7.34  sec.  Notice  that  there  are  four  vortices  visible 
in  the  plot:  two  positive  and  two  reversed.  The  additional 
small  vortices  are  apparently  a  consequence  of  entrainment  by 
the  major  ones.  As  far  as  we  know,  their  existence  has  not 
been  noted  previously. 

When  we  repeat  the  calculation  with  nuclear  initial 
conditions,  several  differences  appear  at  a  very  early  stage. 
The  reflected  shock  propagates  upward  rapidly  through  the 
much  hotter  fireball  and  reverberates  more.  The  maximum  flow 
speeds  (as  opposed  to  sound  speeds)  are  smaller,  a  difference 
which  persists  to  late  times.  Although  the  shock  radius  as  a 
function  of  time  is  essentially  the  same,  the  rarefaction 
wave  moves  faster  as  the  deeper  density  well  gets  filled  in. 

Figures  9(a) -(c)  show  plots  analogous  to  those  of 
Figures  6(g)— (i) 7  by  this  time  it  is  clear  that  much  of  the 
early  difference  in  the  density  profiles  is  washed  out  as 
pressure  begins  to  relax  to  ambient.  (The  pressure  differs 
from  ambient  by  <  5%  everywhere,  so  that  pressure  contour 
plots  are  not  very  informative.)  Note,  however,  the  indenta¬ 
tions  that  appear  on  the  underside  of  the  internal  energy 
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contours  at  the  base  of  the  stem  [Fig.  9(b)].  These  are 
absent  in  the  corresponding  HE  plot.  Fig.  6(g).  We  conjec¬ 
ture  that  they  are  the  signature  of  a  fluid  instability, 

possible  Rayleigh-Taylor .  The  idea  is  that  the  air  sucked  in 
by  the  (forward)  vortex  at  the  bottom  of  the  fireball  is  much 
denser  than  the  fireball  itself.  In  running  into  the  latter 
it  sets  up  the  classic  condition  for  Rayleigh-Taylor  insta¬ 
bility  (direction  of  effective  gravity  and  density  gradient 
are  opposed ) . 

It  is  clear  that  the  major  qualitative  differences 
between  the  HE  and  nuclear  cases  persist  longest  in  the 
velocity  plots.  This  is  not  surprising,  as  the  circulation 
patterns  represented  by  the  vortices  have  essentially 
infinite  lifetimes  in  the  absence  of  viscosity,  we  have  run 
both  nuclear  and  HE  cases  out  to  stabilization  (not  shown 

here)  and  have  shown  that  there  are  qualitative  differences 
in  velocity  plots  to  the  very  end.  At  all  times  t>0  the  peak 
flow  velocity  in  the  HE  Case  exceeds  that  in  the  comparable 
nuclear  result.  This  is  a  reflection  of  the  fact  that  the 
Chapman-Jouguet  solution  at  a  radius  of  10  m  has  a  pressure 
peak  of  52  kbar,  vs  3  kbar  for  the  1-kton  standard  at  the 
same  radius.  The  means  that  the  former  starts  out  with  much 
more  violent  motion,  i.e.,  fluid  velocities  an  order  of 
magnitude  larger.  In  point  of  fact,  the  HE  case  does  not 
closely  resemble  a  point  source.  At  initialization  the 

nuclear  profiles  have  ~  6%  of  the  yield  in  kinetic  energy. 

This  fraction  increases  to  a  maximum  of  ~  15%,  then 

decreases.  In  the  HE  case  the  fraction  is  initially  about 
one-half  and  decreases  monotonically  thereafter  at  about  the 
same  r  ••  .e  as  in  the  nuclear  case. 
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Section  5 
CONCLUSIONS 


We  have  described  numerical  simulations  carried  out 
for  a  600-ton  HE  burst  and  a  1-kton  nuclear  burst,  both  at 
166  ft.  The  code,  gridding,  and  method  of  solution  are  the 
same  in  the  two  calculations.  Although  the  shock  waves  in 
the  two  cases  propagate  at  roughly  the  same  speed  and  break 
away  at  roughly  the  same  time,  and  although  the  pressure 
fields  relax  to  ambient  in  similar  fashion,  we  find  signifi¬ 
cant  differences,  of  which  the  following  appear  to  be  the 
most  important. 

(i)  The  HE  flow  velocities  are  systematically 
larger. 

(ii)  In  the  regular  reflection  region  (underneath 
and  close  to  the  fireball),  the  HE  case 
exhibits  two  overpressure  peaks  at  the 
surface,  rather  th. n  one,  due  to  re-ref lec¬ 
tion  of  the  reflected  shock  from  the  contact 
surface  between  air  and  detonation  products. 

(iii)  For  the  HE  case  the  upjrar  vortex  forms  first, 
followed  by  the  reverse  vortex  near  the  axis 
of  symmetry  and  the  ground.  Adjacent  HE 
products  begin  to  be  entrained  into  a 
positive  vortex  over  a  longer  period  of  time, 
several  seconds.  In  the  nuclear  case  the 
negative  vortex  is  the  dominant  one;  it  is 
larger  then  the  HE,  persists  longer,  and 


contains  larger  velocities  than  the  positive 
vortex  at  comparable  times. 


(iv)  The  HE  flow  establishes  a  pattern  of  four 
vortices,  two  forward  and  two  reversed, 
instead  of  one  of  each. 

(v)  The  stem  of  the  nuclear  fireball  appears  to 
exhibit  a  Rayleigh-Taylor  instability,  absent 
in  the  HE  case. 

Since  the  velocities  on  axis  are  higher  in  the  HE 
case  (the  upper  positive  vortex  is  larger),  fireball  rise  is 
faster  than  in  the  nuclear  case.  The  nuclear  case,  however, 
probably  scours  up  more  dust  because  the  velocities  are 
larger,  and  the  reverse  vortex  is  larger  and  more  persistent. 
It  is  difficult  to  argue  conclusively  on  this  point  because 
so  much  depends  on  terrain,  conditions  in  the  boundary  layer, 
and  other  physical  effects  not  included  in  this  model  (e.g.r 
precursor  heating  and  turbulence).  Further  study  of  the 
tracer  particle  motions  we  have  calculated  is  expected  to  be 
illuminating  in  this  regard. 
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Washington,  DC  20330 
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01CY  Attn  DEX 

01CY  Attn  NTES-S 
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01CY  Attn  0?  981 
01CY  Attn  OP  03EG 


Commander 

Naval  Surface  Weapons  Center 
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01CY  Attn  Drome -WC  (Tech  lib) 


Commander 
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01CY  Attn  library 

Commander 

US  Aray  Foreign  Science  &  Tech  Ctr 
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33 


Director 

BMD  Advanced  technology  Center 

Department  of  the  A rny 
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